#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <cuda.h>

int J_ORDER;
int J_ROW_TEST;
float J_ERROR;
int J_ITE_MAX;

#define N 8 /*Number of blocks*/
#define TPB 1024 /*Threads per block*/


/*
 * Divide os elementos da matriz matrixA pelos valores de sua diagonal.
 * Cada linha eh dividida pelo seu respectivo elemento diagonal. 
 */
__global__ void new_matrixA(float *matrixA, int order) {
        int tid = threadIdx.x + TPB * blockIdx.x;
        // pega o indice do elemento da diagonal
        int diag_idx = (tid / order) * (order + 1);

        if(tid < order * order) { 
                if (tid != diag_idx) {
                        // A[i,j] = A[i,j] / A[i,i]
                        matrixA[tid] /= matrixA[diag_idx];
                }
        }
}

/*
 * Copia para o vetor lineA a linha da matriz matrixA que sera usada para
 * testar o resultado obtido ao final das iteracoes. test_line define a 
 * linha que sera copiada.
 */
__global__ void get_lineA(float *lineA, float *matrixA, int test_line, int order) {
        int tid = threadIdx.x + TPB * blockIdx.x;

        if(tid < order) {
                lineA[tid] = matrixA[tid + order * test_line];
        }
}
        

/*
 * Divide os elementos do vetor dev_vectorB pelas elementos das diagonais da
 * matriz dev_matrixA. Um elemento i qualquer de dev_vectorB eh dividido pelo
 * elemento diagonal da i-esima liunha de dev_matrixA. Em seguida, a funcao
 * faz uma copia de dev_vectorB em dev_vectorX e zera os elementos diagonais
 * de dev_matrixA. 
 */
__global__ void new_vectorB_and_vectorX(float *dev_vectorB, float *dev_matrixA, float *dev_vectorX, float *dev_rB, int order) {
        int tid = threadIdx.x + TPB * blockIdx.x;
        // pega o indice do elemento da diagonal
        int diag_idx = tid * (order + 1);

        if(tid < order) {
                dev_vectorB[tid] /= dev_matrixA[diag_idx];
                dev_vectorX[tid] = dev_vectorB[tid];
                dev_matrixA[diag_idx] = 0.0;
        }
}


/*
 * Faz uma copia do vetor vectorX no vetor old_vectorX, e zera vectorX.
 */
__global__ void updateX(float *vectorX, float *old_vectorX, int order) {
        int tid = threadIdx.x + TPB * blockIdx.x;

        if(tid < order) {
                old_vectorX[tid] = vectorX[tid];
                vectorX[tid] = 0.0;
        }
}


/*
 * Calcula o produto entre a matriz matrixA e o vetor old_vectorX e salva na
 * matriz matrix_prod. Em seguida, soma os elementos de cada linha de matrix_prod
 * e salva no vetor vectorX. A soma da i-esima linha de matrix_prod eh salva na
 * i-esima posicao de vectorX. Por fim, subtrai, posicao por posicao, o valor 
 *  obtido em vectorX dos valores presentes no vetor vectorB, e armazena o 
 * resultado novamente em vectorX. Ao final, o vetor vectorX contem o resultado
 * de uma iteracao completa. 
 */
__global__ void calculateXi(float *vectorX, float *old_vectorX, float *matrixA, float *vectorB, float *matrix_prod, int order) {
        int tid = threadIdx.x + TPB * blockIdx.x; 
        int i;

        if(tid < order * order){
                matrix_prod[tid] = old_vectorX[tid % order] * matrixA[tid];
	}
        // sincroniza todas as threads de um mesmo bloco
        __syncthreads();
        // usa apenas a thread 0 de cada bloco
        if (tid < order) {
        	// faz a reducao da linha bid da matrix_prod na posicao bid de vectorX
                // x[i] = somatorio(A[i][j]*x[j]
                vectorX[tid] = 0.0;
                for (i = 0; i < order; i++) {
                	vectorX[tid] += matrix_prod[tid * order + i];
                }
                // x[i] = B[i] - somatorio(A[i][j]*x[j]
                vectorX[tid] = vectorB[tid] - vectorX[tid];
        }
}


/*
 * Calcula a diferenca entre os vetores vectorX - oldVector e salva 
 * no vetor vectorDiff. Em seguida, calcula a razao entre o maior
 * elemento de vectorDiff e o maior elemento de vectorX e salva no
 * argumento mr. mr vai ser usado para testar se o erro maximo permitido
 * ja foi atingido.
 */ 
__global__ void calculate_mr(float *mr, float *vectorX, float *oldVector, float *vectorDiff, int order){
        int tid = threadIdx.x + TPB * blockIdx.x;
        int i;
        float max_diff, max_vectorX, tmp;

        if(tid < order){
                vectorDiff[tid] = fabs(vectorX[tid] - oldVector[tid]);
	}		
	//sincroniza todas as threads do mesmo bloco 
	__syncthreads();

        // apenas a thread 0 vai executar
        if(tid == 0) {
        	// encontra o maior termo de vectorDiff e o maior modulo de vectorX
                max_diff = vectorDiff[0];
                max_vectorX = fabs(vectorX[0]);
                for (i = 1; i < order; i++) {
                	if(vectorDiff[i] > max_diff) {
                        	max_diff = vectorDiff[i];
                        }
                        tmp = fabs(vectorX[i]);
                        if(tmp > max_vectorX) {
                        	max_vectorX = tmp;
                        }
                }
                *mr = max_diff/max_vectorX;
        }
}


/*
 * Faz o produto entre o vetor lineA e o vetor vectorX e salva
 * no vetor tmp. Em seguida, soma todas as posicoes de tmp e 
 * salva no argumento resp.
 */
__global__ void calculate_res(float *lineA, float *vectorX, float *tmp, float *res, float order) {
        int tid = threadIdx.x + TPB * blockIdx.x;
        int i;

        if (tid < order) {
                tmp[tid] = lineA[tid] * vectorX[tid];
	}
        // sincroniza todas as threads de um mesmo bloco
        __syncthreads();

        // usa apenas a thread 0 
        if (tid == 0) {
        	// faz a reducao do vetor tmp
                // res = somatorio(tmp[i])
                *res = 0.0;
                for (i = 0; i < order; i++) {
                	*res += tmp[i];
                }
        }
}


int main(int argc, char *argv[]){
        float *matrixA;
        float *vectorB;
        float *lineA;
        float *vectorX;
        float *oldVector;
        float *vectorDiff;
        float rB;
        float mr;
        float res;
        int   it = 0;

        float *dev_matrixA;
        float *dev_vectorB;
        float *dev_lineA;
        float *dev_vectorX;
        float *dev_oldVector;
        float *dev_vectorDiff;
        float *dev_rB;
        float *dev_matrix_prod;
        float *dev_mr;
        float *dev_res;

        FILE *file_in = fopen(argv[1], "r");
        int i;
	/* Se falhar ao abrir o arquivo. */
        if( file_in == NULL ){
                printf("FILE ERROR - THE FILE COULDNT BE OPENED\n");
                exit(-1);
        }

	/* Le os parametros do teste */
        rewind( file_in );
        fscanf( file_in, "%d", &J_ORDER); /* ordem da matriz */
        fscanf( file_in, "%d", &J_ROW_TEST); /* linha a ser testada */
        fscanf( file_in, "%f",&J_ERROR); /* erro maximo permitido */
        fscanf( file_in, "%d", &J_ITE_MAX); /* nro maximo de iteracoes a serem executadas */

	/* define o tamanho dos blocos e o tamanho do grid de threads que serao
	   utilizados. Define blocos de TPB threads. O tamanho do grid eh definido
	   de acordo com a ordem O da matriz de entrada. Eh usado N blocos de TPB threads
	   onde N eh o menor nro tal que (N * TPB) > (O * O). */
	dim3 grid(J_ORDER*J_ORDER/TPB + 1);
	dim3 block(TPB, 1);
	/* Aloca memoria para as vairaveis da CPU */
        matrixA = (float*) malloc(J_ORDER*J_ORDER*sizeof(float));
        vectorB = (float*) malloc(J_ORDER*sizeof(float));
        lineA = (float*) malloc(J_ORDER*sizeof(float));
        vectorX = (float*) malloc(J_ORDER*sizeof(float));
        vectorDiff = (float*) malloc(J_ORDER*sizeof(float));
        oldVector = (float*) malloc(J_ORDER*sizeof(float));
	/* Aloca a memoria para as variaveis do device */
        cudaMalloc((void**)&dev_matrixA, J_ORDER*J_ORDER*sizeof(float*));
        cudaMalloc((void**)&dev_vectorX, J_ORDER*sizeof(float));
        cudaMalloc((void**)&dev_vectorB, J_ORDER*sizeof(float));
        cudaMalloc((void**)&dev_oldVector, J_ORDER*sizeof(float));
        cudaMalloc((void**)&dev_lineA, J_ORDER*sizeof(float));
        cudaMalloc((void**)&dev_vectorDiff, J_ORDER*sizeof(float));
        cudaMalloc((void**)&dev_rB,sizeof(float));
        cudaMalloc((void**)&dev_matrix_prod, J_ORDER*J_ORDER*sizeof(float));
        cudaMalloc((void**)&dev_mr,sizeof(float));
        cudaMalloc((void**)&dev_res,sizeof(float));

	/* Le a matriz A */
        for(i=0;i<J_ORDER*J_ORDER; i++){
                fscanf(file_in, "%f", &matrixA[i]);
        }
	/* Le o vetor B */
        for(i=0;i<J_ORDER;i++){
                fscanf(file_in, "%f", &vectorB[i]);
        }
	/* Salva o elemento do vetor B que sera usado no teste ao final das iteracoes */
        rB = vectorB[J_ROW_TEST];

	/* Copia a matriz A e o vetor B para o device */
        cudaMemcpy(dev_matrixA,matrixA,J_ORDER*J_ORDER*sizeof(float),cudaMemcpyHostToDevice);
        cudaMemcpy(dev_vectorB,vectorB,J_ORDER*sizeof(float),cudaMemcpyHostToDevice);

	/* Salva a linha que sera usada no teste ao final das iteracoes */
        get_lineA<<<grid,block>>>(dev_lineA,dev_matrixA,J_ROW_TEST,J_ORDER);
	/* Calcula a matriz A* */
        new_matrixA<<<grid,block>>>(dev_matrixA, J_ORDER);
	/* Calcula o vetor B* e inicializa o vetor X */
        new_vectorB_and_vectorX<<<grid,block>>>(dev_vectorB, dev_matrixA, dev_vectorX, dev_rB, J_ORDER);

        do {
		/* Faz uma copia do vetor X */
                updateX<<<grid,block>>>(dev_vectorX, dev_oldVector, J_ORDER);
		/* Calcula o novo velor do vetor X */
                calculateXi<<<grid,block>>>(dev_vectorX, dev_oldVector, dev_matrixA, dev_vectorB, dev_matrix_prod, J_ORDER);
		/* CAlcula o erro */
                calculate_mr<<<grid,block>>>(dev_mr, dev_vectorX, dev_oldVector, dev_vectorDiff, J_ORDER);
		/* Copia o erro obtido para a CPU */
                cudaMemcpy(&mr,dev_mr,sizeof(float),cudaMemcpyDeviceToHost);
        } while (mr > J_ERROR && ++it < J_ITE_MAX);
	/* Calcula o valor da equacao presente na linha de teste de acordo com o resultado obtido
	   nas iteracoes acima. */
        calculate_res<<<grid,block>>>(dev_lineA, dev_vectorX, dev_oldVector, dev_res, J_ORDER);
	/* Copia o resultado da equacao da linha de teste para a CPU */
        cudaMemcpy(&res,dev_res,sizeof(float),cudaMemcpyDeviceToHost);

        printf("---------------------------------------------------------\n");
        printf("Iterations: %d\n", it);
        printf("RowTest: %d => [%f] =? %f\n", J_ROW_TEST, res, rB);
        printf("---------------------------------------------------------\n");

	/* Libera a memoria do device*/
        cudaFree(dev_matrixA);
        cudaFree(dev_vectorB);
        cudaFree(dev_vectorX);
        cudaFree(dev_vectorDiff);
        cudaFree(dev_oldVector);
        cudaFree(dev_rB);
        cudaFree(dev_lineA);
        cudaFree(dev_matrix_prod);
        cudaFree(dev_mr);
        cudaFree(dev_res);
	/* Libera a memoria da CPU */
        free(matrixA);
        free(vectorX);
        free(vectorB);
        free(vectorDiff);
        free(oldVector);
        free(lineA);

        return 0;
}